Check out the ReadMe page on Github for a description of JulES. This is a demo of the prototype implementation of JulES in a stepwise manner. Here are some of the most important steps:
using DataFrames, Plots, Statistics, JSON, Distributed
plotlyjs();
The WebIO Jupyter extension was not detected. See the WebIO Jupyter integration documentation for more information.
The problem is simulated on 31 2.20 GHz processors running in parallel. TODO: Test on faster processors.
const numcores = 31
if nprocs() < numcores
addprocs(numcores - nprocs())
end
@show nprocs();
nprocs() = 31
@everywhere include(joinpath(dirname(pwd()),"jgrc/TuLiPa/src/TuLiPa.jl"));
@everywhere include(joinpath(dirname(pwd()),"jgrc/JulES/src/JulES.jl"));
datayear = 2025
scenarioyearstart = 1981
scenarioyearstop = 2011
totalscen = 30 # scenarios to consider uncertainty for
# Standard time for market clearing - perfect information so simple time type
tnormal = FixedDataTwoTime(getisoyearstart(datayear),getisoyearstart(scenarioyearstart))
# Phasein settings
phaseinoffsetdays = 2 # also simulation step length
phaseinoffset = Millisecond(Day(phaseinoffsetdays)) # phase in straight away from second stage scenarios
phaseindelta = Millisecond(Week(5)) # Phase in the second stage scenario over 5 weeks
phaseinsteps = 5 # Phase in second stage scenario in 5 steps
# Make scenario times for all uncertainty scenarios. List of tuples with tnormal, tphasein and scenarionumber
totalscentimes = []
for scen in 1:totalscen
scentnormal = FixedDataTwoTime(getisoyearstart(datayear),getisoyearstart(scenarioyearstart + scen - 1))
scentphasein = PhaseinTwoTime(getisoyearstart(datayear),getisoyearstart(scenarioyearstart), getisoyearstart(scenarioyearstart + scen - 1), phaseinoffset, phaseindelta, phaseinsteps);
push!(totalscentimes, (scentnormal, scentphasein, scen))
end
This is a simplified version of the dataset NVE uses for its long-term power market analyses (for example https://www.nve.no/energi/analyser-og-statistikk/langsiktig-kraftmarkedsanalyse/). The dataset consist of:
We cannot publish the full dataset, but many of the profiles for wind, solar, hydro and demand can be downloaded from https://www.nve.no/energi/analyser-og-statistikk/vaerdatasett-for-kraftsystemmodellene/.
sti_themadata = "data_fra_thema"
themastructure = json_to_elements(sti_themadata, "dataset_thema.json")
themaseries = json_to_elements(sti_themadata, "tidsserier_thema.json")
elements = vcat(themaseries,themastructure)
addscenariotimeperiod!(elements, scenarioyearstart, scenarioyearstop);
# Set horizons for price prognosis models
# Long
longhorizonduration = Millisecond(Week(5*52))
longhydroperiodduration = Millisecond(Week(6))
longrhsdata = StaticRHSAHData("Power", datayear, scenarioyearstart, scenarioyearstop)
longmethod = KMeansAHMethod()
longclusters = 4
longunitduration = Millisecond(Hour(6))
# Medium
medhorizonduration = Millisecond(Week(54))
medhydroperiodduration = Millisecond(Day(7)); @assert medhorizonduration.value % longhydroperiodduration.value == 0
medrhsdata = DynamicRHSAHData("Power")
medmethod = KMeansAHMethod()
medclusters = 4
medunitduration = Millisecond(Hour(4))
# Short
shorthorizonduration = Millisecond(Week(1))
shorthydroperiodduration = Millisecond(Day(1)); @assert medhorizonduration.value % shorthorizonduration.value == 0
shortpowerparts = 12
# Preallocate storage for problems and results on different cores. Use package DistributedArrays
# Distribute scenarios
allscenarios = distribute(totalscentimes)
# Problems are built, updated, solved, and stored on a specific core. Moving a problem between cores is expensive, so we want it to only exist on one core.
longprobs = distribute([HiGHS_Prob() for i in 1:length(allscenarios)], allscenarios)
medprobs = distribute([HiGHS_Prob() for i in 1:length(allscenarios)], allscenarios)
shortprobs = distribute([HiGHS_Prob() for i in 1:length(allscenarios)], allscenarios)
# Results are moved between cores. These are much smaller than longprobs/medprobs/shortprobs and are inexpensive to move between cores.
medprices = distribute([Dict() for i in 1:length(allscenarios)], allscenarios)
shortprices = distribute([Dict() for i in 1:length(allscenarios)], allscenarios)
medendvaluesobjs = distribute([EndValues() for i in 1:length(allscenarios)], allscenarios)
nonstoragestates = distribute([Dict{StateVariableInfo, Float64}() for i in 1:length(allscenarios)], allscenarios)
# Organise inputs and outputs
probs = (longprobs, medprobs, shortprobs)
longinput = (longhorizonduration, longhydroperiodduration, longrhsdata, longmethod, longclusters, longunitduration, "long")
medinput = (medhorizonduration, medhydroperiodduration, medrhsdata, medmethod, medclusters, medunitduration, "long")
shortinput = (shorthorizonduration, shorthydroperiodduration, shortpowerparts, "short")
allinput = (numcores, elements, allscenarios, phaseinoffset)
output = (medprices, shortprices, medendvaluesobjs, nonstoragestates)
# Initialize price prognosis models and run for first time step. Run scenarios in parallell
@time pl_prognosis_init!(probs, allinput, longinput, medinput, shortinput, output)
86.040697 seconds (5.53 M allocations: 252.401 MiB, 1.22% compilation time)
Task (done) @0x00000000744e3650
scenario = 1
index = collect(medprices[scenario]["steprange"])
prices = medprices[scenario]["matrix"]
labels = [name for name in medprices[scenario]["names"]]
p = plot(index,prices,label=reshape(labels, 1, length(labels)),legend=:outertopright)
for scenario in 2:5 # length(allscenarios)
prices = medprices[scenario]["matrix"]
labels = [name for name in medprices[scenario]["names"]]
plot!(p,index,prices,label=reshape(labels, 1, length(labels)),legend=:outertopright)
end
display(p)